N2O valence (orbs 6 - 9) MF-PADs with polarization¶

18/03/24 PH

N2O with polarization dependence (this notebook).

N2O basic LF results: https://phockett.github.io/ePSdata/N2O-preliminary/N2O_orbs6-9_LF_preliminary_220224-tidy.html

TODO:

  • Rerun in updated env and better plots.
  • Skip slurm files! Currently these are also named .out by default - should check header.
  • NormInt checks not working correctly here? Change in output format?
  • Automatic reordering for dir read? By orb number.

For methods: https://epsproc.readthedocs.io/en/dev/demos/ePSproc_class_demo_161020.html

In [ ]:
# Quick hack to override default HTML template
# NOT required in some JLab versions.
# https://stackoverflow.com/a/63777508
# https://stackoverflow.com/questions/21971449/how-do-i-increase-the-cell-width-of-the-jupyter-ipython-notebook-in-my-browser
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))

Setup¶

In [1]:
import sys
import os
from pathlib import Path
import numpy as np
# import epsproc as ep
import xarray as xr

import matplotlib.pyplot as plt

from datetime import datetime as dt
timeString = dt.now()

import epsproc as ep

# Plotters
from epsproc.plot import hvPlotters

# Multijob class dev code
from epsproc.classes.multiJob import ePSmultiJob

hvPlotters.setPlotters(width = 700, snsStyle='whitegrid')

OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

* sparse not found, sparse matrix forms not available. 
* natsort not found, some sorting functions not available. 

* Setting plotter defaults with epsproc.basicPlotters.setPlotters(). Run directly to modify, or change options in local env.

* Set Holoviews with bokeh.
* pyevtk not found, VTK export not available. 

* Set Holoviews with bokeh.

Load data¶

In [2]:
import warnings
# warnings.filterwarnings('once')   # Skip repeated numpy deprecation warnings in current build (xr15 env)
# warnings.filterwarnings('ignore')   # Skip repeated numpy deprecation warnings in current build (xr15 env)
In [3]:
# # Scan for subdirs, based on existing routine in getFiles()

# fileBase = Path('/home/paul/ePS/OCS/OCS_survey')  # Data dir on Stimpy

# fileBase = Path('/home/paul/fock-mount/globalhome/eps/N2O/N2O_valence') # Data dir on Jake (fock-mount)
# fileBase = Path('/home/paul/ePS/N2O/N2O_valence') # Data dir on Jake (local copy)
fileBase = Path('/home/jovyan/jake-home/ePS/N2O/N2O_valence') # Data dir on Jake (from ePSprc container)
In [4]:
# TODO: fix orb label here, currently relies on (different) fixed format

data = ePSmultiJob(fileBase, verbose = 0)

data.scanFiles()
data.jobsSummary()

*** Warning: Missing records, expected 32, found 24.
*** Warning: Found 8 blank sets of matrix elements, symmetries ['A2']
*** Warning: Missing records, expected 32, found 24.
*** Warning: Found 8 blank sets of matrix elements, symmetries ['A2']
*** Warning: Missing records, expected 32, found 24.
*** Warning: Found 8 blank sets of matrix elements, symmetries ['A2']
*** Warning: Missing records, expected 32, found 24.
*** Warning: Found 8 blank sets of matrix elements, symmetries ['A2']
*** Warning: Missing records, expected 32, found 24.
*** Warning: Found 8 blank sets of matrix elements, symmetries ['A2']
*** Warning: Missing records, expected 32, found 24.
*** Warning: Found 8 blank sets of matrix elements, symmetries ['A2']
*** Warning: Missing records, expected 32, found 24.
*** Warning: Found 8 blank sets of matrix elements, symmetries ['A2']
*** Warning: Missing records, expected 32, found 24.
*** Warning: Found 8 blank sets of matrix elements, symmetries ['A2']
Found 4 directories, with 16 files.

*** Job orb9_P details
Key: orb9_P
Dir /home/jovyan/jake-home/ePS/N2O/N2O_valence/orb9_P, 4 file(s).
{   'batch': 'ePS N2O, batch N2O_valence, orbital orb9_P',
    'event': 'orb 9 (P/CAv) ionization, batch N2O_valence, None.',
    'orbE': -13.499568716634647,
    'orbLabel': 'orb9_P',
    'symmetryLabel': 'P/CAv'}

*** Job orb6_S details
Key: orb6_S
Dir /home/jovyan/jake-home/ePS/N2O/N2O_valence/orb6_S, 4 file(s).
{   'batch': 'ePS N2O, batch N2O_valence, orbital orb6_S',
    'event': 'orb 6 (S/CAv) ionization, batch N2O_valence, None.',
    'orbE': -22.60177741591763,
    'orbLabel': 'orb6_S',
    'symmetryLabel': 'S/CAv'}

*** Job orb7_P details
Key: orb7_P
Dir /home/jovyan/jake-home/ePS/N2O/N2O_valence/orb7_P, 4 file(s).
{   'batch': 'ePS N2O, batch N2O_valence, orbital orb7_P',
    'event': 'orb 7 (P/CAv) ionization, batch N2O_valence, None.',
    'orbE': -21.11331458826209,
    'orbLabel': 'orb7_P',
    'symmetryLabel': 'P/CAv'}

*** Job orb8_S details
Key: orb8_S
Dir /home/jovyan/jake-home/ePS/N2O/N2O_valence/orb8_S, 4 file(s).
{   'batch': 'ePS N2O, batch N2O_valence, orbital orb8_S',
    'event': 'orb 8 (S/CAv) ionization, batch N2O_valence, None.',
    'orbE': -19.00715329282262,
    'orbLabel': 'orb8_S',
    'symmetryLabel': 'S/CAv'}

*** Job stacked details
Key: stacked
No 'job' info set for self.data[stacked].

System properties¶

Note orbital numbering in table below:

  • Orb is Gamess file output orbital numbering, but energy but not grouped by degeneracy.
  • OrbGrp is grouped numbering by degeneracy, also used by ePolyScat, and will be used for labels etc. below.

BUT - file names are 0-indexed (oops), so give OrbGrp-1 here. Sorry.

In [5]:
data.molSummary()

*** Molecular structure

*** Molecular orbital list (from ePS output file)
EH = Energy (Hartrees), E = Energy (eV), NOrbGrp, OrbGrp, GrpDegen = degeneracies and corresponding orbital numbering by group in ePS, NormInt = single centre expansion convergence (should be ~1.0).
props Sym SymOrb EH Occ E NOrbGrp OrbGrp GrpDegen NormInt
orb
1 S 0.0 -20.6612 2.0 -562.219894 1.0 1.0 1.0 0.975090
2 S 0.0 -15.8503 2.0 -431.308635 1.0 2.0 1.0 0.975090
3 S 0.0 -15.7018 2.0 -427.267745 1.0 3.0 1.0 0.975090
4 S 0.0 -1.6338 2.0 -44.457963 1.0 4.0 1.0 0.975090
5 S 0.0 -1.4519 2.0 -39.508212 1.0 5.0 1.0 0.975090
6 S 0.0 -0.8306 2.0 -22.601777 1.0 6.0 1.0 0.975090
7 P 0.0 -0.7759 2.0 -21.113315 1.0 7.0 2.0 0.975090
8 P 0.0 -0.7759 2.0 -21.113315 2.0 7.0 2.0 0.999998
9 S 0.0 -0.6985 2.0 -19.007153 1.0 8.0 1.0 0.975090
10 P 0.0 -0.4961 2.0 -13.499569 1.0 9.0 2.0 0.975090
11 P 0.0 -0.4961 2.0 -13.499569 2.0 9.0 2.0 0.999998 
*** Warning: some orbital convergences outside single-center expansion convergence tolerance (0.01):
[[ 7.  2.]
 [ 8.  2.]
 [10.  2.]
 [11.  2.]]

Compute MF-$\beta_{L,M}$ as function of ellipticity¶

Define field(s)¶

Using EfiledPol class, set various ellipticities and fields $(E_x,E_y)$. For MF calcs, compute for all polarization geometries corresponding to fields aligned along $(x,y,z)$ molecular axes.

  • Frame definitions docs
  • Methods: EfieldPol class demo
In [6]:
from epsproc.efield.epol import EfieldPol
import numpy as np

# Mulitple ellipticities and rotations
# For fields, set with format = [amplitude,azimuth (0,pi), ellipticity (-pi/4,pi/4)]
states = 10

# angles = np.linspace(0, np.pi, states)  # For range of angles
angles = np.zeros(states)  # no rotation

maxEll = 0.5  # Set max ellipticity to use (0 - 1)
ellipticities = np.linspace(0,maxEll*np.pi/4, states)   # Set for 0 - max ellipticity
ell = np.c_[np.ones(states),angles,ellipticities]
labels = (ellipticities/(np.pi/4) * 100).round(2)  # Set labels at %age ellipticity

# Set fields
Eell_multi = EfieldPol(ell=ell, labels = labels)

# Plot fields
Eell_multi.plot(figsize=(6,6), draw_arrow=False)
plt.legend(labels, loc='upper right')

Set field from ell.
[[1.         0.         0.        ]
 [1.         0.         0.04363323]
 [1.         0.         0.08726646]
 [1.         0.         0.13089969]
 [1.         0.         0.17453293]
 [1.         0.         0.21816616]
 [1.         0.         0.26179939]
 [1.         0.         0.30543262]
 [1.         0.         0.34906585]
 [1.         0.         0.39269908]]
Out[6]:
<matplotlib.legend.Legend at 0x7fc8d30405e0>

Compute MF results¶

In [8]:
# Compute MF results
data.MFBLM(EfieldPol = Eell_multi, thres = 1e-1)

Set parameters to `self.epDict` and `self.epXR`.
Set pol state data to self.YLM and self.YLMrot, and orientations to self.RX.
Set orientations to self.epDict['RX'].
Set geomCalc.EPR() results to `self.EPRX`.

Plot all results vs. Eke¶

Note pol slider is ellipticity %age as set above (0 - 50%).

In [9]:
# Plot BLM parameters
# Note this may be quite slow for large results sets
data.BLMplot(dataType='MFBLM', thres = 1e-2, backend='hv', width=700,  ylim=(-1.5, 2.5), xlim=(0,15))

BLMplot set data and plots to self.plots['BLMplot']
Out[9]:
True

Plot all results vs. ellipticity¶

Note pol axis is ellipticity %age as set above (0 - 50%).

In [10]:
# Plot BLM parameters - VERSUS POL STATE
# Note this may be quite slow for large results sets
# To speed up pass subselection options to plotter, or replot (see next cell).
data.BLMplot(dataType='MFBLM', xDim='pol', Erange=[0,10], thres = 1e-2, backend='hv', width=700,  ylim=(-1.5, 2.5)) #, xlim=(0,15))

BLMplot set data and plots to self.plots['BLMplot']
Out[10]:
True

Plot all results vs. ellipticity and layout orbitals¶

Note pol axis is ellipticity %age as set above (0 - 50%).

In [18]:
# Example with post-plotter subselection for faster plotting...

# Create HV object, but don't plot
data.BLMplot(dataType='MFBLM', xDim='pol', thres = 1e-2, backend='hv', width=1000,  ylim=(-1.5, 2.5), 
             renderPlot=False)

# Replot from Holoviews dataset with selectors
# Note this may need hv or hvPlotters to be loaded
from epsproc.plot import hvPlotters
xDim = 'pol'
# data.plots['BLMplot']['hvDS'].select(Eke=[1.1,10.1,20.1]).select(Orb='orb5').to(hvPlotters.hv.Curve, kdims=xDim).overlay(['l', 'm'])
data.plots['BLMplot']['hvDS'].select(Eke=[1.1,10.1,20.1]).to(hvPlotters.hv.Curve, kdims=xDim).overlay(['l', 'm']).opts(width=1000,  ylim=(-1.5, 2.5)).layout('Orb').cols(1)

BLMplot set data and plots to self.plots['BLMplot']
Out[18]:
In [14]:
# Quick data IO - https://epsproc.readthedocs.io/en/3d-afpad-dev/dataStructures/ePSproc_dataStructures_IO_demo_280622.html#Pickle
# see Pemtk for wrappers

# import pickle
# # Save an Xarray
# with open(Path(fileBase, 'N2O_MF-polStateScans_180324.pickle'), 'wb') as handle:
#     pickle.dump(data.data, handle, protocol=pickle.HIGHEST_PROTOCOL)

# ep.IO.writeXarray
In [ ]:

Versions¶

In [19]:
!hostname

8fdb9101787f
In [20]:
import scooby
scooby.Report(additional=['epsproc', 'xarray', 'jupyter'])
Out[20]:
Mon Mar 18 16:07:33 2024 EDT
OS Linux CPU(s) 64 Machine x86_64 Architecture 64bit
RAM 62.8 GiB Environment Jupyter File system btrfs
Python 3.10.11 | packaged by conda-forge | (main, May 10 2023, 18:58:44) [GCC 11.3.0]
epsproc 1.3.2-dev xarray 2022.3.0 jupyter Version unknown numpy 1.23.5
scipy 1.10.1 IPython 8.13.2 matplotlib 3.5.3 scooby 0.7.2
In [21]:
# Check current Git commit for local ePSproc version
!git -C {Path(ep.__file__).parent} branch
!git -C {Path(ep.__file__).parent} log --format="%H" -n 1

* 3d-AFPAD-dev
565448918fb2dfcb8b5860b0174f5711d85edd2c
In [22]:
# Check current remote commits
!git ls-remote --heads https://github.com/phockett/ePSproc

565448918fb2dfcb8b5860b0174f5711d85edd2c	refs/heads/3d-AFPAD-dev
7e4270370d66df44c334675ac487c87d702408da	refs/heads/dev
1c0b8fd409648f07c85f4f20628b5ea7627e0c4e	refs/heads/master
69cd89ce5bc0ad6d465a4bd8df6fba15d3fd1aee	refs/heads/numba-tests
ea30878c842f09d525fbf39fa269fa2302a13b57	refs/heads/revert-9-master
baf0be0c962e8ab3c3df57c8f70f0e939f99cbd7	refs/heads/testDev
In [ ]: